Stress-induced patterns in ion-irradiated Silicon: a model based on 

anisotropic plastic flow 

Scott A. Norris 
Department of Mathematics 
Southern Methodist University, Dallas, TX 7527^\ 
(Dated: July 23, 2012) 

Abstract 

We present a model for the effect of stress on thin amorphous films that develop atop ion-irradiated 
silicon, based on the mechanism of ion-induced anisotropic plastic flow. Using only parameters 
directly measured or known to high accuracy, the model exhibits remarkably good agreement with 
the wavelengths of experimentally-observed patterns, and agrees qualitatively with limited data 
on ripple propagation speed. The predictions of the model are discussed in the context of other 
mechanisms recently theorized to explain the wavelengths, including extensive comparison with an 
alternate model of stress. 
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I. INTRODUCTION 



Among the many nanoscale patterns that have been observed on ion-irradiated surfaces 
[TJ |2], the discovery of hexagonal arrays of high- aspect ratio dots on Ar + -irradiated GaSb 
|3j has sparked a flurry of experimental and theoretical study into spontaneous pattern 
formation as a potential route to "bottom-up" fabrication of nanoscale devices. A growing 
body of evidence increasingly suggests that highly-ordered structures occur only for targets 
containing more than one material [IH6]. Although early observations of dot formation 
on pure Silicon [7H9] sparked a series of proposed explanations for the ordered structures 
1 10, 11 j, none turned out to be viable [12, 1 3 J : ultimately, dots disappeared when impurities 
and geometric artifacts were carefully removed [H], and reappeared upon their systematic 
re-introduction pQ, [T5HT7] . As a result, much recent attention has focused on theories of 
irradiated binary materials [6j [T81-23J . 

Despite these results on the formation of ordered structures, the study of monatomic tar- 
gets remains important because of commonalities between binary and monatomic systems, 
and the relative simplicity of the latter. In particular, because it is readily amenable to 
molecular dynamics simulation, and its near-surface region is amorphous under ion bom- 
bardment, noble-gas ion irradiation of silicon has been an important system for comparison 
between experiment, theory, and simulation of pure materials. This allows rapid develop- 
ment and testing of theories on the basic physical processes of ion irradiation, which occur on 
two time scales. On timescales of individual ion impacts, spanning ~ 10~ 9 sec and described 
as the prompt regime, there is erosion of some target atoms away from the target [21H26] . 
and redistribution of others to new locations [TT| 127]. On much longer timescales associated 
with kinetic relaxation, spanning around ~ 10 2 sec and denoted the gradual regime, there 
is the possibility of surface diffusion [26] , stress buildup [28], and viscous flow [29]. The un- 
derstanding of these basic processes, which remain important in binary materials, continues 
to evolve in important ways, bringing closer the goal of a theory developed enough to make 
accurate predictions on irradiation-induced structure formation. 

In this paper, we investigate the effect on morphology evolution of stress buildup and 
relaxation in the thin film of material affected by ion irradiation. We choose a different 
approach than a recent series of papers on this topic [301432] . which cast the role of stress 
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in terms of classical fluid dynamics results on gravity-driven flows |33j ; the relationship 



between that model and our own will be discussed in detail in Section IV Here, instead, we 
draw analogy with the large literature on ion-induced stress for irradiation at energies in the 
MeV range, where incoming ions are primarily slowed by electronic stopping [34TI39] . There, 
incoming ions induce an anisotropic deformation in the material, leading to stress buildup. 
More recently, this effect has also been observed during irradiation in the keV range [40j, 
and the mathematical form of the model has been successfully used to describe a number 
of phenomena at even lower energies where nuclear stopping is dominant |4T| |4*2"] . Based on 
these results, we here apply the model formally to pattern formation at low energies. We find 
that our model is able to predict the observed wavelength to remarkable accuracy without 
the need to estimate any unknown parameter values. We also find it makes predictions on 
the ripple propagation speed that are consistent with experimental observations. 



II. MODEL 



We consider the morphological evolution of the top layer of a monatomic target, irradi- 
ated by noble gas ions at an incidence angle 9. The ions penetrate a certain average distance 
into the solid before initiating a collision cascade of atomic displacements. Over many im- 
pacts, the accumulated damage induced by these cascades amorphizes a thin film of material 
atop the target, with a crystalline /amorphous boundary at the bottom the film, and a free 
boundary at the top of the film. We construct a co-ordinate system in which the (planar) 
target is perpendicular to the z-axis, with z = h(x, y) tracking the top free boundary, and 
z = g (x, y) tracking the bottom boundary. We will then consider the effect of stress on the 
stability of this configuration. In this paper, we will focus exclusively on the effect of stress, 
and therefore neglect the prompt regime of sputtered atoms [24, 26J, or those redistributed 
to new locations [TT1 1271 133] . 



A. Constitutive Law 



Basic Material Properties. It has been established for some time that the incoming ions 
impart an effective fluidity to the amorphous film that is many orders of magnitude larger 
than that of bulk amorphous silicon [29] - So any treatment of the film must include viscous 
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effects. Proposals in past years have suggested the additional consideration of elastic effects 
[HI HI], and a constitutive law including all such effects is |39j 

E = — T + 4^T+ — — (trT)I + E B (1) 

2r] 2G 9B Dt v ' B K ' 

Here the first three terms contributing to the linear rate-of-strain tensor 

E = \ (Vv + Vv r ) (2) 



constitute a standard linear model for a Maxwell fluid [44J, where T is the stress tensor, 
r) is the viscosity, G is the shear modulus, and B is the bulk modulus. The final term on 
the right, E_b, is a contribution to this rate of strain due the ion beam. However, it has 
been argued elsewhere that, in a variety of relevant regimes, the contribution of elasticity 
should be negligible [T3| 142] . Hence, we take the limit of a purely viscous {G — > oo) and 
incompressible {B — > oo) fluid. In the incompressible limit, however, a hydrostatic pressure 
term must also be included, leading to the simpler constitutive law we will consider for the 
remainder of this work, 

T = —pi + 277 (E - E B ) , (3) 

which is nearly equivalent to a Newtonian fluid, except for the addition of the constant 
contribution —2r]EiB to the stress. 

Effect of Stress. We now turn to the modeling of stress. Recent work by Castro, Cuerno, 
and co-workers has appealed to analogy with familiar concepts in fluid dynamics by assuming 
that the microscopic mechanism stress generation can be coarse-grained into an effective body 



force [30| |3T] . We will discuss that approach in Section IV, but will here take a different 
approach, attempting to connect more directly with the microscopic process. At that scale, 
MD simulation has shown that each impact significantly redistributes the target silicon atoms 
to new locations |43"| H5J 0B], gradually increasing the magnitude of a compressive stress to 
a saturated state |H], which is also observed experimentally [48J. Each impact thus induces 
a direct deformation of the material, suggesting that the effect of the beam be incorporated 
directly into the constitutive relationship between stress and strain, in a way that depends 
linearly on the total fluence. 

A model with exactly these properties has already been developed to describe anisotropic 
plastic flow during high-energy ion irradiation in the regime of electronic stopping, where 



rapid thermal cycling due to ion impacts leads to a stress-free strain rate of the form j3] 



E B = fAD (0) . 



(4) 



Here, / is the ion flux, A is a measure of the magnitude of strain induced per ion, and D 
describes the (anisotropic and angle-dependent) shape of that strain. The mechanism of be- 
hind Eqn.Q is of course not directly applicable in the nuclear stopping regime. However, the 
phenomenon anisotropic plastic flow has been observed even in the nuclear stopping regime 
|40| . and the mathematical form of the governing equations has been applied successfully to 
describe various phenomena at low energies |1TI WI\ . In this spirit, we have argued elsewhere 
that, for = 0, simple symmetry arguments lead almost immediately to the form 



D(0) 



1 
1 
0-2 



(5) 



regardless of the underlying mechanism. From here, it is natural to add angle-dependence 
via the matrix 

cos — sin (0) 
R(0) = 10 • (6) 
sin(0) cos(0) 

describing rotation about the y-axis, which leads to 



D(0) = R(-0)D(O)R(0) 



cos (20) - | 



| sin (20) 




sin (20) -§ cos (20) - \ 



Combining (|3]), Q, we obtain the final form of the constitutive relation, 



(7) 



pi + 2rj E - fAD 



with D (0) given by Eqn. ([?]). 
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B. Governing Equations 



The governing equations are obtained following standard continuum analysis in the limit 
of a highly viscous film. In the bulk, conservation of mass in the incompressible limit gives 

V-v = 0, (8) 

the incompressibility condition on the velocity field v. In addition, conservation of momen- 
tum takes the form 

P (W +V ' VV ) = V ' T + f > ( 9 ) 
where the left hand side contains the density p and the material acceleration, and the right 
hand side lists forces V-T due to internal stresses, and f due to long-range body forces. From 
here, because of the high viscosity of the ion-irradiated film, accelerations in the left-hand 
side are small and can be neglected, and we do not consider the effect of any body forces. 
Hence, we arrive at Stokes equations: 

V ■ T = -Vp + r]V\ = 0. (10) 

It remains to apply boundary conditions. At the amorphous/crystalline boundary z = 
g (x, y) we have the combined "no-penetration" and "no-slip" conditions for a viscous fluid; 

v = 0. (11) 

Meanwhile, allowing for the effect of surface tension on a viscous fluid, the stress balance at 
the free boundary z = h(x,y) can be written 

T ■ n = -7/cn (12) 

while the kinematic condition relating bulk velocity v to the normal surface velocity at the 
interface, Vj, is 

Vj — v • n (13) 
We see that the effect of the beam, through the angle-dependent tensor D(0), appears 



mathematically only at the free boundary in Equation (12), where it alters the conditions 
necessary for stress balance. 
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III. ANALYSIS 



A. Steady Solution 



We first look for a steady state (d/dt — > 0) consisting of a flat, uniform film occupying 
the space between z = go and z = ho, with pressure po and velocity vo- For convenience we 
choose go = 0. Symmetry considerations greatly simplify the calculations - the assumption 
of uniformity implies translational symmetry in x and y, which limits the steady pressure 
and velocity field to be functions of z at most: p (z) and v (-2). Furthermore, although 
off-normal ion incidence from the negative x- direction may break reflection symmetry in x, 
we enforce reflection symmetry in y, which prohibits a velocity component in the ^/-direction. 
These considerations immediately limit the form of the steady velocity field Vo to 

v„(z) = Mz),0, Wo{z)) T . (14) 



Comparing Eqn. (14) to the incompressibility condition (pj) and the no-penetration condition 



(11), we first see that 



w (z) = 0. 



(15) 



Next, the ^-component of the Stokes equations (10) show the steady pressure to be a constant 

p (z) = a, (16) 



while the x-component of ( 10 ) restricts uq (z) to at most the linear form 

uo (z) = b + cz. (17) 
Finally, application of the remaining boundary conditions gives us a, b, and c. The no-slip 



condition (11), applied at z = go = 0, requires 6 = 0, while the stress balance condition (12) 



applied at z = ho, yields a = rjfA [3 cos (29) + 1] and c = 3 / A sin (29). In summary, the 
steady solution is thus 

Po = vfA [3 cos (29) + 1] 
v = (3 f Asm (29) z, 0, 0) T 

cos (29) 

cos 2 (9) 





-6 V fA 
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consisting of a uniform pressure and downbeam shear that both depend on the irradiation 
angle-of-incidence 9. The pressure is maximum at 9 = 0° and monotonically decreases 
at higher angles, whereas the shear flow is maximal at 9 = 45° and decays to zero for 
9 = {0°, 90°}. We note in particular the form of the steady stress when 9 = 0; a symmetric, 
biaxial compressive stress with no vertical component. 



B. Linear Stability 



We now perform a linear stability analysis - e.g., we consider the evolution of the system 
when the top boundaries is subjected to an infinitesimal, sinusoidal perturbation in the x- 
and y- directions. This leads to solutions that can be expressed in terms of normal modes 
having the form 
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7) i(k 1 x+k2y)+crt 



(19) 



Here we do not assume that the bottom boundary remains flat. Instead, as the top boundary 
changes shape, the zone of material affected by collision cascades changes with it. As a 
simplest approximation to this behavior, based on the observed weak dependency of film 
thickness on angle [49J, we choose 

9i = hi 

indicating that the bottom boundary is simply a vertical translation of the top boundary, 



as depicted in Figure III C With this definition, the forms (19) are then inserted inserted 



into the governing equations ( 10 )-( 13 ), and terms are kept only to leading order in the 



vanishingly-small parameter e. Solution proceeds in three steps, as shown in the Appendix. 



First, the linearized version of the incompressibility condition (\8h and Stokes equations (10) 



are solved to give a general solution for the pressure and velocity fields in the bulk. Second, 



the linearized boundary conditions (11) and (12) are applied to determine the resulting 



integration constants. Last, having uniquely determined the pressure and velocity fields, the 



application of the kinematic condition (13) leads directly to the dispersion relation a (ki, k 2 ) 



governing the evolution of each mode. We find 
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a(k 1 ,k 2 



-QfA 



cos (29) (hoh) 2 + cos 2 (9) (h k 2 ) 
1 + 2Q 2 + cosh (2Q) 



-3/Asin (26)i{hoh) 



,2 cosh (Q) [Q 2 + sinh 2 (Q)l / s 

3/4 sin (20) 2 (/loifei) { r , ^\ J - cosh (Q) 



1 + 2Q 2 + cosh (2Q) 



(20) 



7 Q (sinh (2Q) - 2Q) 



2rjh 1 + 2Q 2 + cosh (2Q) 
where the dimensionless wavevector Q = h ^kf + fcf. Here the first line corresponds to the 
effect of anisotropic plastic flow operating on the bulk film, the second line corresponds to 
the effect (under said flow) of the non-planar bottom boundary gi (x, y) = hi (x, y), and the 
third line corresponds to the well-studied (isotropic) effect of surface leveling under surface 
tension (see Ref. |50|). 



A commonly-employed simplification of Eqn.(20) occurs in the long- wavelength limit Q 



1 of Eqn.(20), in which wavelength is much longer than the film depth. In this limit the 



leading-order contributions of each line in Eqn.(20) reduce to 



a w -3fA [cos (29) (Mo) 2 + cos 2 (9) (k 2 h ) 2 ] 
-pAsm(29)t(k 1 h )Q 2 



(21) 
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3r)h, 



-Q 4 



This approximation is useful by nature of being much simpler than Eqn.(20). However 



despite being commonly employed, this assumption is less commonly verified against exper- 
imental data. For instance, in the experiments of Madi [91 [51], widely used as a point of 
comparison with theory, the film thickness h pa 3nm, while wavenumbers are as large as 
k p^ 2n/ (20 nm), giving Q Pd 1. Hence, while we include Eqn.(21) for reference and use in 



qualitative discussion, numerical predictions will be based on the full Eqn.(20). 



C. Discussion of Results 



Stability. The real terms 



r (k) = Re [a (k)} 



from Eqn.(21) describe the growth rate of perturbations to the flat steady state in terms of 



their wavenumber k. The value of k that maximizes r (k) - denoted k* - is called the most 



unstable mode; if a (k*) < 0, then all modes decay and the steady solution is stable , whereas 
if a(k*) > 0, then at least some modes near k* grow, and the steady solution is unstable. 



From the full dispersion relation (20), we see that a positive band of unstable wavenumbers 



exists if 9 > 45°. In the longwave approximation (21 ), the most unstable wavelength is given 

by 



A* « 2tt, 



2jh 



2vr, 



47/1 



(22) 



9fArj cos (20) V 3 |T (0) | cos (29) 
where the denominator inside the radical is identified as a multiple of the steady in-plane 



stress at normal incidence, as given by Eqn.( 18). This is convenient, because it allows replac- 



ing the estimation of the infrequently-measured ion-enhanced viscosity 7], and the stress-per- 
ion parameter A, with an experimental observation of the steady stress at normal incidence. 
Hence, we can compare existing measurements of the steady stress and the wavelength for 
the same material under identical irradiation conditions. It has already been shown by Madi 
that at 250 eV, h « 3nm and |T (0)| ~ 1.5 GPa 1481 . while in the same chamber at the 



same energy, wavelengths given by Figure III C were observed |2J dH EI]- With the addi 



tional value 7 = 1.36 J/m [52| 153] . we can compare the stress measurements to the observed 
wavelengths. We find that the full dispersion relation predicts a remarkably good fit to the 
data. Unsurprisingly given the size of Q, the longwave approximation does not do nearly as 
well. 

Ripple Velocity. The imaginary terms 

00 (k) = -Im [a (k)} 



from Eqn.(21) (note the conventional minus sign) describe the propagation rate of perturba- 
tions to the flat steady state because, from the original linear ansatz, 



exp (at + ikx) = exp (rt + ik(x — —t) 



(23) 



Equation (23) reveals that, in a stationary "lab" frame of reference, individual ripples of 



wavelength k propagate with a speed 



V, 



phase 



u(k) 



10 




Predicted Wavelength 



— full solution 

- - LW approx. 




50 



60 
angle 



70 



80 



90 



Figure 1. (color online) Wavelength as predicted by numerically maximizing the full dispersion 
relation ( |20[ ) over Q (solid line), and under the longwave approximation ( |21[ ) (dashed line). Given 
the fixed parameter values indicated in the text, the full dispersion relation predicts the observed 
wavelengths of parallel-mode ripples (circles) to surprising accuracy, although it does not predict the 
rotation of the ripples to perpendicular-mode (square) at high angles. Following the discussion the 
text, we observe that the wavelengths predicted by the longwave approximation differ significantly 
from those predicted by the full dispersion relation. 



usually called the phase velocity. Comparison with (21) reveals that, for our problem 



. ,->ro,h(Q) [Q 2 + sinh 2 (Q)] 

W = SfAh* sin (26) < 1 + 2g2 + cosh(2g) " ^ h (Q) 
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fAh sin (29) Q 5 



(24) 



Although limited experimental data on ripple velocities is available, what reports exist indi- 
cate that - contrary to early seminal theory concerned only with erosion [26J - ripples always 



propagate in the direction of the ion beam |54"t 155] . Both our longwave result ( 24 ) , and also 



the velocity resulting from the full dispersion relation (20), are uniformly positive for all 
values of 6 or Q, consistent with these observations. 
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Figure 2. Schematic of the system, illustrating the intuitive understanding of ripple growth and 
propagation. 



Intuitive Explanation of Ripple Growth and Propagation. Both the instability and trans- 
lation mechanisms can be understood intuitively in terms of simple differences in the steady 
shear flow and pressure across a ripple. The instability occurs when the net shear flow on 
the uphill side exceeds the net flow on the downhill side. We recall that the flow rate has a 
maximum when the local angle of incidence is 45°, and decreases monotonically for as the 
local angle continues to increase. So, when 9 > 45, then uphill slopes (with angle closer 
to 45°) have greater flow than downhill slopes (with higher angle), leading to mass accu- 
mulation at the top of the ripple. Similarly, ripple translation occurs because the steady 
pressure is a monotonically decreasing function of incidence angle. Hence, for any angle, the 
pressure in the uphill side exceeds that in the downhill side, leading to flow from uphill to 
downhill slopes, as indicated in Figure [Til C[ We note that this ripple propagation mechanism 



is almost entirely driven by the nonplanar nature of the bottom boundary; if the bottom 
boundary were instead flat, the only contribution to the phase velocity would come from the 



last term of the first line of (20), which indicates simple translation along with the shearing 
film. However, when the translated bottom boundary is included, it exactly cancels this 
simple shear in the long- wavelength limit, and the resulting ripple propagation is entirely 
due to differences in pressure across a ripple. 
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Relationship to thermal stress models. The contribution (24) to the ripple velocity, and 
indeed several aspects of the analysis used to attain it, bears a notable resemblance to the 
work of Alkemade on ripple propulsion There, the author proposes a mechanism of 

ripple propagation that also is caused by variations in the stress, but in which the stress 
itself is caused by thermal expansion. Briefly, larger relative depth of the uphill slopes 
causes higher sustained temperatures there, leading to greater thermal expansion and thus 
greater thermal stress. The problem is that at the scale of the ripple size (~ 20 — 50 nm), 
thermal diffusion is so rapid that no meaningful temperature differential can be sustained. 
Our model can be thought of as an alternative to this mechanism, where the source of stress 
is purely mechanical, requiring no thermal gradients, and variation appears naturally as a 
result of the different slopes inherent to a ripple structure. 

IV. RELATIONSHIP TO EFFECTIVE BODY FORCE MODELS 

We will now attempt to give a comprehensive comparison of our theory with another 
recent model of stress that appeals to analogies with fluid dynamics [301432] by rendering 
the effect of stress as an "effective body force" (EBF) which acts throughout the amorphous 
film. The bulk of the the EBF model is developed in Ref.|31j. There, the effect of the beam 
is posited to lie within the constitutive relation via 

Tij = -pSij + fi (dM + djVt) + , (25) 

which is effectively identical to our Eqn.Q, and similar to work dating at least to Ref.[39j. 
However, the authors never explicitly define the tensor describing the effect of the beam, 
instead proposing only that its divergence takes the form 

V • T s = b = f E V (6 - 7) (26) 

where the authors "define b = V • T s as a body force acting in the bulk of the fluid layer." 
Although the incoming ions undoubtedly exert a force on the film as they are slowed down by 



the film, that force can be shown to be vanishingly small [49J, and Eqn.(26) is not suggested 
to represent any actual physical force. Instead, f# is proposed to contain "the coarse-grained 
information about the effect of the residual stress created in the target ... [with] dimensions 
of a gradient of stress," and {6 — 7) to encode dependence upon the local angle of incidence 
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6 — 7- This is an interesting effort to encompass the many complexities of the nonlocal ion 
irradiation process within a familiar, intuitive form, and seems to achieve good agreement 
with experiment. It therefore demands reasons for the proposal of an alternative model. 



A. Differences with EBF: Fundamental Model 

We here present three concerns with the EBF model that are avoided by the present 
approach. 

Surface-Dependent Body Force. The first, primary concern with the EBF model given 



by Eqn. ( 26 ) is its highly unusual suggestion that a "body force" could exist that acts through- 
out the film, but which contains an explicit dependence on the configuration of some nearby 
patch of surface. Traditionally, true body forces that can act directly on bulk elements - also 
known more descriptively as "long-range forces" - are completely independent of nearby sur- 
faces (since the latter are presumably much nearer than the source of the long-range force). 
Any effects due to surfaces are instead rigorously segregated into boundary conditions, to 
ensure that the transmission of surface physics into the bulk is mediated by the constitutive 



properties of the material itself. A formulation like (26) blurs this traditional distinction, 
suggesting the existence of a mechanism which depends on the surface, but can nevertheless 
directly inform buried bulk parcels without any consideration of the material in between. 

Functional Form of^. A second concern is that, to obtain a good fit with experimental 
data, the EBF model assumes that the function \l/ take the form \1/ = cos (8 — 7); i.e., that 
the body force is proportional to the flux density at the surface. This would be perfectly 
reasonable if the force b acted only at the surface - i.e., if it were a surface traction - 
because it reflects the geometric dilution of the beam as the orientation of the surface varies 
relative to the beam orientation. However, the force b is instead proposed to act as a long- 
range force throughout the bulk of the film, and as just discussed, true long-range forces 
are typically approximately constant in any given location. Hence, the correct form for the 
angle dependence of \l/ is much less obvious than it first appears, and the plausibility of the 
assumption ^ = cos (0 — 7), upon which good agreement with experiment rests, becomes 
less certain. 
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Location of Stress. A final, minor concern is the inclusion of the undefined tensor T s 
into the stress T. True body forces are traditionally not included within the stress tensor 
itself, but rather constitute a separate item in the equations of conservation of momentum 
(J9j). The distinction is subtle: the resulting equations for the velocity field are the same in 
either case, but once velocities have been calculated, the resulting stress field ^ depends 
qualitatively on the location of the effects of the beam. So within the EBF model, if the effect 



of the beam is genuinely supposed to appear within Eqn.(25), then T s must be rigorously 
defined because, as a part of the full stress T, it must play a role in the balance of forces at 
the free interface z = h(x,y). However, if the effect of the beam is truly to be treated as a 



body force (as seems to be the intent), then it should not appear in Eqn.(25). 

Comparison with the present model. In contrast to these dilemmas inherent in the EBF 
model as currently proposed, the model presented here maintains the traditional location and 
form of terms, and makes fewer assumptions to achieve good agreement with experiment. In 
particular, the ion beam acts within the bulk constitutive law Q as a true tensor describing 
stress-free anisotropic plastic flow. Acting as it does throughout the entire body of the film, it 
is completely independent of the surface configuration, as a traditional long-range mechanism 
should be. And, as a true tensor, it is fully defined via Eqn.([7|) and included in the force 



balance (12) at the free boundary. Finally, by starting at a level closer to the microscopic 
process, our model based on the mechanism of anisotropic plastic flow requires fewer, and 
weaker, assumptions. In particular, the correct angular dependence of the wavelength on the 
irradiation angle arises very naturally merely from the simple rotation of the strain tensor 
with the beam. 



B. Differences with EBF: Theoretical Predictions 

In addition to fundamental differences at the level of modeling, the present theory also 
differs from the EBF theory in terms of experimental predictions. 

Steady State. The two models for stress produce different steady states. We first observe 
that, in contrast to the constant pressure and linear flow profile found here, Ref.|3T] reports 
a depth-dependent pressure and quadratic flow profile. This results in qualitatively different 
scalings for the net flow through the film: d 2 vs. d 3 , which could be significant given the 
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small size of d. Next, combining Eqs.(4),(15-17) from Ref.|31|. one can show that the EBF 
model results in a steady stress throughout the film of 



TW = f E * {0) z 



(27) 



cos (9) - sin (0) 
— sin (9) cos (9) 

where the second component T s remains undefined within the EBF theory, and should 
potentially be omitted, as discussed above. Focusing on first part associated with viscous 



stress, we observe two differences with our Eqn.(18). First, the magnitude of the viscous 



stress increases as the distance from the free surface increases, reaching a maximum at the 
crystalline/amorphous interface z = —d. However, considering that this stress is due to 
incoming ions which are slowed by the film, one would expect the effect of the beam to 



at most remain constant with film depth (as assumed here in Eqn.(18)), or, perhaps more 



realistically, to decay with film depth. Second, Eqn.(27) has a different shape than our 



Eqn.(18), with different angular dependence and principal stress directions, which may be 



detectable using wafer curvature measurements. In particular, when dotted with the upward- 



pointing unit vector k, Eqn.(27) indicates the presence of stress along the vertical direction; 



however, unless a body force is truly present, the flat free boundary should allow all stress 



in the vertical direction to relax to zero (as seen here in Eqn.(18)). 

Wavelengths. A noteworthy success of the EBF theory is the apparently remarkable 



agreement of its predictions on wavelength with experimental data; our Figure ( III C ) is 



nearly identical to Figure 3 of Ref . |31J . However, it turns out that the remarkable agreement 
with experiment depicted in that figure were based on incorrect data. Taking the geometric 
mean of MD estimates of 1.62 GPa (56], and preliminary experimental data of 200 MPa 
|51j . the authors estimate the magnitude of the stress to be 569 MPa. This value then 
leads to an estimate of |fg| = 0.424 kg/(nms) 2 , producing the good fit with data. However, 
the lower, experimental value of the stress was later found to be in error, and a more 
correct value for Silicon irradiated by Argon at 250 eV has more recently been revised to 
1.4 GPa |48j. Although the precise means of calculating |fe| are not specified in Ref. |31|. if 
we assume that |fg| oc |To| we may estimate a value associated with the revised stress of 
\¥e\ ~ 0.424 x 1400 -j- 569 = 1.04kg/(nms) 2 , which does not fit the data on wavelengths as 
well. Put another way, the EBF model seems to over-estimate the effect of stress on growth 
rate by a factor of about 2-3 relative to the current theory, which provides better agreement 
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with experiment when correct values of the steady stress are used. 

Ripple Velocities. The theory of solid flow developed around the EBF model in Ref. 
predicts a ripple velocity of the form 

d 2 f E sm(29) / 3_ 2 



^ppie = ' [I ~ -Q 2 ) (28) 

(although the derivation appears to contain a few minor errors, such as the use of the group 



velocity ^ instead of the phase velocity % ). At first glance, Eqn.(28) appears to have 



the same sin (29) angle dependence as our result 24; however, it also bears a few important 



differences. First, within the parenthesis, Eqn.(28) contains a constant term, and because 
it is typically assumed that Q C 1 (though see above), this term should be expected to 



dominate the ripple velocity for large wavelengths. In contrast, our Eqn.(24) contains no 
constant term which, because the critical value of Q itself depends on 9, should produce 
a different functional form of the ripple velocity in terms of the angle. This difference 
should be especially noticeable near the bifurcation at 45°, where wavelengths become large: 
our model predicts stationary ripples as Q — > 0, whereas the EBF model predicts moving 



ripples. Second, Eqn.(28) depends on the ion-enhanced viscosity of the material, but only 
weakly or not at all on the flux (see Ref.[32j for a discussion of parameter dependencies 



of Je)- In contrast, our Eqn.(24) does not depend on the viscosity, and depends linearly 
on the flux, consistent with the hypothesis that each impact induces a fixed amount of 
anisotropic deformation in the material. These differences should be testable experimentally 
with measurements of ripple velocity as functions of various experimental parameters. 



V. OPEN QUESTIONS 

Ripple Rotation. One of the fundamental successes of the original Bradley-Harper the- 
ory of erosion was its ability to predict the oft-observed 90-degree rotation of ripples at high 
angles of incidence |26j. Despite recent studies suggesting that this erosive effect is in fact 
overwhelmed at most angles by the effect of redistributed atoms |43l [57] , no alternative to 
Bradley-Harper has been able to definitively predict ripple rotation at high angles, and our 
theory is no exception. It is likely that either erosive effects are still somehow dominant 
at grazing angle, or that some other effect particular to grazing incidence induces the rip- 
ple rotation. (We point out that Johnson et al. recently predicted such a rotation based 
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on molecular dynamics measurements [58]. However, that prediction rests upon apparent 
measurements of a negative lateral momentum transfer at grazing angles; we have instead 
observed lateral momentum transfer to be positive at all angles [43J. Comparative study of 
the methods used to obtain and filter these data would be highly desirable.) 

Relationship to Mass Redistribution. A notable observation absent from the discussion 
surrounding the EBF model is the observation that stress is actually the second distinct phys- 
ical mechanism shown to produce reasonably good agreement with experimentally-observed 
ripple wavelengths. It has been shown previously that the leading order treatment of 're- 
distributed' target atoms - those displaced from their original location but not sputtered 



- produces second derivatives in the dispersion relation (21) of exactly the same form as 
those reported here |11[ [27] , and molecular dynamics simulations suggest that these terms 
are of just the right size to produce reasonably accurate predictions of the observed ripple 
wavelengths |13]. Because each effect, alone, seems able to predict the correct wavelengths, 
the combination of both effects would be expected to be wrong. 

One possible resolution to this dilemma appeals to the fact that stress, ultimately, is 
generated by target atom displacements. Hence, morphology evolution due to both target 
atom redistribution and stress are ultimately driven solely by single ion impacts, just on 
different timescales. So although they are properly distinct within a timescale-separated 
framework, one could imagine the possibility of a generalized multi-scale framework which 
incorporates both the displacements and stresses induced by single ion impacts, as measured 
by molecular dynamics. In the process of analyzing the contribution to stress (see, e.g., 
Refs.|47 [ l56 | l59|). the tensor D is something that could potentially be measured as a function 
of angle, whereas the effective body force, until it is connected directly with a microscopic 
mechanism, will likely be more difficult to obtain. 

A second possible explanation appeals to uncertainty surrounding the exact value of the 
ion-enhanced viscosity r], which has only been estimated using molecular dynamics. Whereas 
the wavelength predictions in Ref.[43j due to mass redistribution relied directly on 77, within 
the present approach it was possible to avoid estimating rj by instead measuring the steady 
stress I Tq I, which can be directly observed experimentally. Hence, while the accuracy of 
the redistribution model is vulnerable to changes in estimates of r), the accuracy of the 
anisotropic plastic flow model is not. For this reason, the latter currently appears to this 
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author to be the more robust. So despite its technical problems, the authors of the EBF 
model may ultimately be proved correct in their hypothesis that "Solid flow drives surface 
nano-patterning by ion-beam irradiation" [32J. 



VI. CONCLUSIONS 



In the pursuit of better understanding the role of stress in ion-irradiated films, we have 
analyzed a model of stress based on the mechanism of anisotropic plastic flow. This approach 
represents an alternative to a recent "effective body force" (EBF) model of stress based 
loosely on traditional fluid mechanics results. Based on distinct microscopic mechanism, 
the model avoids several technical problems in the EBF approach, and achieves remarkable 
experimental agreement on ripple wavelengths using more accurate experimental parameters 
and fewer underlying assumptions. In addition, it makes more intuitive predictions on the 
steady stress, and alternate, testable predictions on the ripple velocity. Although challenges 
remain, notably a failure to predict ripple rotation at grazing incidence angles, and a need to 
distinguish stress effects from those of redistributed atoms, this closer link to physical origins, 
and expanded predictive power for Argon-irradiated Silicon, make the model a potentially 
valuable stepping stone toward the goal of a predictive theory with the quantitative accuracy 
necessary for industrial relevance. 
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Appendix A: Solution of the Linearized Equations 



In this Appendix we present the details of the calculation by which the dispersion relation 



(20) was obtained. 



General solution from the bulk Equations. In the linear ansatz (19), the (already-linear) 
Stokes equations become 

- ik x p (z) +r}[- (fc? + k l) u + u"] = 
-ik 2 p (z) + (k\ + kf)v + v"] = 
-p (z) + n[- (k\ + k 2 2 )w + w") = 
ik\U + ik 2 v + w' = 

The first step in the solution is to eliminate u,v,w in the following manner: 



1. differentiate Eqn.(A4), to get an expression for w" 



2. insert the resulting expression for w" into Eqn.(A3) 



3. differentiate the result to get an expression for p" (z) 



4. in the resulting expression, eliminate w' using (A4) 



(Al) 
(A2) 
(A3) 
(A4) 



5. finally, eliminate u and v using (Al) and (A2) 
The resulting homogeneous equation for p (z) is 

P" ~ (*? + kl) p = 

with solutions 

p(z) = A cosh (Rz) + B sinh (Rz) , (A5) 

where R = \Jk\ + k\. 

From here, we simply insert the solution (A5) back into (Al)-(A3) to get inhomogeneous 
equations for u, v, w: these have solutions 



V 



w 



ik\ 



C cosh (Rz) + D sinh (Rz) H -z \B cosh (Rz) + A sinh (Rz 



2t]R 

ik 2 



E cosh (Rz) + E sinh (Rz) + — -z [B cosh (Rz) + A sinh (R 
G cosh (Rz) + H sinh (Rz) + —z [A cosh (Rz) + B sinh (Ri 



(A6) 
(AT) 
(A8) 
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At this point, however, we have too many constants because we first differentiated Eqn.(A4) 



raising its order. Re-considering Eqn.( |A4| in its original form, we insert ( A6 )-( A8 ) and collect 
like terms in sinh (Rz) and cosh (Rz) to obtain 



^D + ^F + ^-B + G 
R R 2r]R 



iki 



C 



iko 



-E 



-A + H 




0. 



Ml iZ) 



v 1 [z 



Wi (z 



C cosh (Rz) + D sinh (Rz) - 
E cosh (Rz) + F sinh (Rz) - 
G cosh (Rz) + H sinh (Rz) - 



R 

ik 2 
~~R 



R R 2rjR 

From here, we choose to replace A and B with the other constants to obtain the following 
expressions for the pressure and velocity: 

Pl ( z ) = -2r] [(RH + ikxC + ik 2 E) cosh (Rz) + (RG + ik x D + ik 2 F) sinh (Rz)] 

ik x (RG + ik x D + ik 2 F) z cosh (Rz) 
+ (RH + ikxC + ik 2 E) z sinh (Rz) 

(RG + ik x D + ik 2 F) z cosh (Rz) 
+ (RH + ihC + ik 2 E) z sinh (Rz) 

(RG + ik x D + ik 2 F) z sinh (Rz) 
\- (RH + ikiC + ik 2 E) z cosh (Rz) 

where R = ^/k\ + kj, and C, D, E, F, G, H are integration constants. 

Integration constants from Boundary Conditions. To find the integration constants, we 



(A9) 



apply the linearized versions of the three no penetration/slip conditions (11) at the amor- 



phous/crystalline interface z — hi (x,y), and the three stress balance conditions (12) at the 



free interface z = ho + hi (x,y). At the bottom of the film, z = gi (x,y) = hi (x,y), the 



no-slip boundary condition (11) linearizes to the related condition 



at z = 0, giving 



vi (x, y, 0) + ^ (x, y,0)-h 1 = 

C = -3fAsin(26) hi (x,y) 
E = 
G = 



(MO) 



(All) 



At the top of the film, z — ho + hi (x,y), we apply the stress condition (12); again, this 
linearizes to a related condition 

T • n x + Ti • n = -7«in 
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at z = ho, which forms a matrix equation for the remaining constants{D, F, H}: 



RC + ^(C + 2QS) ^(C + 2QS) -2ihQC 
^(C + 2QS) RC + §(C + 2QS) -2ik 2 QC 
-2ik x QC -2%k 2 QC 2R (C - QS) 

where Q = h R, C = cosh (Q), S = sinh (Q), and 



+ 3r}fAht sin (26) 





D 




a 




F 




13 




H 




7 



(A12) 



a 









%k\ cos (26) 




= -7 





+ &rjfAh! 


ik 2 cos 2 (6) 


7 




R 2 








RS + | (S + 2QC) 

^{S + 2QC) 

—2ik\Q sinh (Q) 

(A13) 



The three column vectors on the right hand side of Eqn.(Al3) represent, respectively, the 
effects of surface tension, the beam stress acting in the bulk, and the effect under that stress 
of the non-planar bottom boundary. 



The solution to Eqn.(Al2) can be found using Cramer's Rule or a computer algebra 
system, yielding the result 

D = i [2aR 2 (C 2 - QSC) + 2k 2 (ak 2 - pk x ) (C 2 + QSC + 2Q 2 ) + 2i^k x RQC 2 \ 



A 
1 

A 
1 



[2f3R 2 (C 2 - QSC) - 2h (ak 2 - f3h) (C 2 + QSC + 2Q 2 ) + 2i^k 2 RQC 2 } , (A14) 



H = — [2iR (ah + f3k 2 ) QC 2 + 2jR 2 (C 2 + QSC)] 



where a, (3, 7 remain as given in Eqn.(Al3) and 



A = 2R 3 cosh (Q) [1 + 2Q 2 + cosh (2Q)] 



is the determinant of the matrix on the left hand side of Eqn.(Al2) 



Dispersion Relation from the Kinematic Condition. Having uniquely determined the 



pressure and velocity fields, it remains to apply the kinematic condition (13), which linearizes 
to 



a (ki, k 2 ) = wi (h ) - u (h ) - v (h ) ^" 1 



(A15) 



dx u v " y dy 

From the general solution (A9) in the main text, and the solutions (All) and ( A14[ ) for the 
integration constants, we can evaluate W\ (ho), which leads directly to the dispersion relation 



(20) in the main text. 
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